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ABSTRACT 

We introduce a new method for treating Comptonization in computational fluid dynam¬ 
ics. By construction, this method conserves the number of photons. Whereas the traditional 
“blackbody Comptonization” approach assumes that the radiation is locally a perfect black- 
body and therefore uses a single parameter, the radiation temperature, to describe the ra¬ 
diation, the new “photon-conserving Comptonization” approach treats the photon gas as a 
Bose-Einstein fluid and keeps track of both the radiation temperature and the photon number 
density. We have implemented photon-conserving Comptonization in the general relativistic 
radiation magnetohydrodynamical code KORAL and we describe its impact on simulations of 
mildly super-critical black hole accretion disks. We And that blackbody Comptonization un¬ 
derestimates the gas and radiation temperature by up to a factor of two compared to photon- 
conserving Comptonization. This discrepancy could be serious when computing spectra. The 
photon-conserving simulation indicates that the spectral color correction factor of the escap¬ 
ing radiation in the funnel region of the disk could be as large as 5. 

Key words: accretion, accretion discs - black hole physics - relativistic processes - methods: 
numerical 


1 INTRODUCTION 

Accretion on black holes (BHs) is an important energy source that 
shapes the evolution of the Universe. Accreting supermassive BHs 
in the centers of galaxies expel a significant fraction of their ac¬ 
creted rest mass energy in the form of kinetic, thermal and radiative 
feedback, which affects the evolution of their host galaxies (e.g., 
|Fabian|207^|King & Pounds |2015| l. At the same time, stellar mass 
BHs accreting in compact binaries, although not powerful enough 
to influence their hosts, are often the brightest X-ray sources in their 
galaxies. 

For long our understanding of accretion was based on ana¬ 
lytical models, most importantly the thin disk model of IShakur*^ 
|& Sunya^ ( |1973| ) and [Novikov & Thome| ( |1973] ). Although such 
models provide important information on the dynamical properties 
of the accretion disk and predict its radiative luminosity, they pro¬ 
vide only a rough idea of the emitted electromagnetic spectrum. 
To calculate the spectrum properly, including the effect of spec¬ 
tral hardening due to Comptonization in the strongly scattering 
medium, one has to solve the full radiative transfer problem for the 
disk atmosphere. Valuable work has been done in this area by [Davis j 
[et al.[ p005[ ), who computed the vertical structure of disk atmo¬ 
spheres with a self-consistent treatment of radiative transfer. Their 
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work made use of the code TLUSTY ( [Hubeny & Lanz[19^, with 
modifications to treat relativistic disks as described in Hubeny &[ 
[Hubeny [ ( [1997[ ) and [Hubeny et al.[ ( [2001[ ). Spectra obtained with this 
approach have been successfully used for modelling disk emission 
in the thermal-dominant (high/soft) state of X-ray binaries, and in 
particular for estimating BH spins via the continuum fitting method 
(e.g., [McClintock et al.[201l| l. 


The thin disk model can be applied only to optically thick, 
geometrically thin accretion flows. Once the accretion rate exceeds 
the Eddington rate, the various assumptions behind the thin disk 
model (e.g., radiative efficiency, Keplerian rotation) break down. 
Super-critical accretion flows have been modeled semi-analytically 
( [Abramowicz et al.[1988[[S^dowski[20TT] ), but these models do not 
consistently account for mass outflows. In addition, although super¬ 
critical disks are known to be geometrically thick, their spectra have 
in the past been calculated using the plane-parallel approximation 
(see [Straub et al.[20TT] ), which is a gross simplifcation. 


The proper, although computationally expensive, approach 
to modeling accretion disks is to simulate the accretion flow in 
multiple dimension, including angular momentum transport via 
magnetically-driven turbulence { [Balbus & Hawley|[1998[ ), and to 
solve in parallel the multi-dimensional radiative transfer problem. 
Numerical codes to carry out such computations have become 
available only recently. Because of computational cost, often the 
radiative transfer in such simulations is simplified, and only in 
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the postprocessing stage is the full, frequency dependent radia¬ 
tive tranport computed. The first magnetohydrodynamical numeri¬ 
cal simulations of optically thick accretion disks, including radia¬ 
tion, were performed by |Ohsuga et al.| ( [2009| ) and |Ohsuga & Mi-] 
|neshige|(2011| ). These authors used a Newtonian approximation for 
gravity and modeled radiation under the flux limited diffusion ap¬ 
proximation. In recent years, other groups have developed methods 
for general relativistic radiation-MHD simulations, including bet¬ 
ter radiation closures ( [Sqdowski et al.|20l3l [McKinney et al.|2013[ 
[Fragile et al.|2014| l. Newtonian simulations with advanced radiative 
transport methods have also been carried out ( | Jiang et al.|2014aj l. 

Super-critical accretion flows around stellar-mass BHs are 
very optically thick to electron scattering. They have temperatures 
of order a keV or more, so their absorption opacity is subdomi¬ 
nant to scattering. As a result one expects significant Comptoniza- 
tion and spectral hardening of the escaping radiation. So far, how¬ 
ever, simulations of super-critical disks have either ignored Comp- 
tonization altogether ( jS^dowski et al.|20T4} [Jiang et al.|2014bj l or 
have treated it in a crude manner, which we call “blackbody Comp- 
tonization”, in which the radiation is assumed to be locally a perfect 
blackbody ( [Kawashima et al.[2009[ [Sqdowski et al.[2015a[ t. 

It might be hoped that it is sufficient to carry out simula¬ 
tions with the above “blackbody Comptonization,” followed by a 
more detailed proper Comptonization calculation in a later post¬ 
processing stage. However, as we show in this paper, the thermal 
state of the gas depends strongly on the particular implementation 
of Comptonization. The blackbody Comptonization approximation 
should therefore be treated with caution. 

In this work we introduce an improved model of Comptoniza¬ 
tion in which we avoid the assumption that the local radiation at 
each point has a perfect blackbody spectrum. In the new method, 
we independently evolve both the radiation energy density and 
the photon number density, thereby relaxing the previous assump¬ 
tion that the two quantities are related precisely by the Planck 
function. To assess what effect this improvement in the model¬ 
ing of Comptonization has, we have performed three simulations 
of a super-critical accretion disk, one with our previous “black¬ 
body Comptonization,” one with the improved “photon-conserving 
Comptonization” described in this paper, and one with no Comp¬ 
tonization at all. We compare the properties of the three simulations 
and draw conclusions. 

The paper is organized as follows. In Section we introduce 
the blackbody and photon-conserving versions of Comptoniza¬ 
tion. In Section we give implementation details of the photon- 
conserving approach. In Section|^we present the three simulations 
and analyze the results, and in Section|^we briefly discuss the im¬ 
plications of the results. 


Here is the radiation four-force density acting on the gas. In or¬ 
thonormal coordinates in the fluid (or gas) frame (denoted by hats), 
and in the absence of Comptonization, we have 

& = K^{E-AnB\ (4) 

G‘ = {K, + K,,)pF‘. (5) 

Here is the Rosseland mean absorption opacity (cross-sections 
per unit mass), which we adopt for simplicity instead of the proper 
absorption and Planck mean opacities, is the scattering opacity, 
E is the radiative energy density in the fluid frame, B is the intensity 
of blackbody radiation for gas of temperature Tg, 


where a is the radiation constant, and is the radiative flux three- 
vector. 

Equation 0 describes the rate of change of the fii^d energy 
density as a result of energy gain through absorption, /CapE, and en¬ 
ergy loss through emission, /Cap(4;r5). The simplest approximation 
for modeling radiation is to treat it as a local blackbody at each 
point, so that the energy distribution is fully described by the lo¬ 
cal radiation temperature T^. In this approximation, the radiation 
energy density is given by 

E = aT*. (7) 

Correspondingly, we rewrite equation 0 as 

& = K^a(T* - T*). ( 8 ) 

Equation ([^ shows that gas gains energy at a rate proportional to 
and loses energy proportional to Tg. As a result, the two tem¬ 
peratures are pushed towards each other, i.e., the system is driven 
towards thermal equilibrium. Note that k^s does not appear in the 
energy equation. Eor the present discussion, which ignores Comp¬ 
tonization, we have pure Thomson scattering, which only redirects 
photons but does not change photon energies. 

Equation 0 describes the rate of change of the fluid momen¬ 
tum density. The gas acquires the momentum of each photon that 
it either absorbs or scatters. This explains the presence of both 
and /Ces in this equation. However, radiation that is emitted or scat¬ 
tered by the gas is symmetric in the fluid frame. Therefore, there 
is no corresponding momentum loss, hence no counter-balancing 
term with a negative sign. Note that the fluid gains momentum den¬ 
sity in a direction parallel to the radiation flux three-vector F\ and 
the radiation loses a corresponding amount of momentum density. 
Hence the system is driven towards a state in which there is no rela¬ 
tive motion between the fluid and radiation frames, i.e., no radiation 
flux in the fluid frame. 


2 COMPTONIZATION 

2.1 Basic equations of GRRMHD 

Within the framework of general relativity, the interaction between 
gas and radiation is described through the following general rel¬ 
ativistic radiation magnetohydrodynamics (GRRMHD) conserva¬ 
tion equations for gas density, p, and the str ess-energy tensors, Ty 
and Ry, of the magnetized gas and radiation |S^dowski et al.l2013 
[2014] ), 


= 0, 

(1) 

(TX = 

(2) 

= -Gv. 

(3) 


2.2 Thermal Comptonization: Blackbody Approximation 

The main effect of Comptonization is that scattering causes not 
just momentum transfer between the radiation and gas, but also 
energy transfer. A soft photon of energy eq which scatters off a 
thermal electron with temperature on average gains an energy 
(Ae) given by the following expressions in the non-relativistic and 
ultra-relativistic limits, respectively ( [Rybicki & Lightman[1979[ ), 

<Ae> = (9) 

\m^c^ j 

(4kT 

(Ae) = - ^ 6o, 4kT^ » m^c^. (10) 

\m^c^ I 
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For a general temperature, the result is given in equation (2.43) 
in |Pozdnyakov et al.| ( |198!3| . A good fitting function (maximum 
fractional error 1.2%) is the following: 


(A6> = 6o 



1 + 3.683 



kT, 



( 11 ) 

The above expressions are valid so long as the photon is soft, 
i.e., the radiation temperature is much less than the gas temperature. 
When the two temperatures are equal, thermodynamics guarantees 
that there is no energy transfer between gas and radiation. Similarly, 
when the radiation temperature is larger than the gas temperature, 
we expect energy to fiow from the radiation to the gas. To allow for 
these effects, we modify the expression for in equation to the 
following: 


& = paK,(T^-Tp-K,,pE 


, ( 12 ) 

where we have replaced Te by Tg. The negative sign in the Compton 
term is because gas cools when Tg> T^. The cooling is proportional 
to the radiation energy density, E, and to the number of scatterings 
per unit time, ck^sP- The final two factors in square parentheses cor¬ 
respond to the approximate correction for relativistic temperatures 
given in equation |TT] ). Apart from these factors, equ ation (T^ is 
identical to the prescription used by |Kawashima et ar] ( |2009j L 

As far as the momentum equation is concerned, we assume 
that the Compton-scattered radiation is symmetric in the fluid frame 
and carries no net momentum (a fairly good approximation in the 
soft photon limit). Under this approximation, we keep equation ([^ 
unchanged. 

Thus, in this (simplest) version of thermal Comptonization, 
which we call “blackbody Comptonization,” we solve equations ^ 
and using equations and 1^, where Tg and are given by 
equations 0 and Q. respectively. The main weakness of this ap¬ 
proach is that it does not conserve photon number during scattering. 
Instead it assumes that the radiation is a perfect blackbody and uses 
equation 0 to obtain the radiation temperature. 


1 -h 3.683 


kT, 


+ 4 


kT^ 


1 -h 


kT^ 


^k(Tg - T,) 

MgC^ 


2.3 Photon-conserving Comptonization 

In the next approximation, which we call “photon-conserving 
Comptonization,” we give up the assumption that the spectral shape 
is a perfect blackbody, but instead assume that the photons have a 
Bose-Einstein (BE) distribution in the fluid frame, described by two 
parameters: radiation temperature Tr and radiation chemical poten¬ 
tial jd. Correspondingly, we keep track of not only the radiation 
energy and momentum density, E, F\ but also the photon number 
density, Tf. This generalization allows us to include the effect of ra¬ 
diation dilution and spectral hardening. 

Let us define the dimensionless photon energy, x = hv/kT^, 
and the dimensionless chemical potential, ^ = filkT^. The angle- 
integrated BE distribution f(x)dx of photons at temperature is 
given by 

f{x)dx = C{kT,fe~^ x e ’ ^ (^ 3 ) 

ex — e k 

The radiation energy density and photon number density are then 

E = C{kTrfe-dE(4), h{0= - -7, (14) 

Jo e* - e-^ 


r^dx 

- - 

0 


(15) 


In the Planck limit, i.e., when ^ = 0, we have h - 6.493940, 
In - 2 .404114, while in the Wien limit, i.e., when ^ » 1, we have 
Ie — In — 2. ^ ^ 

Noting that ^ Efn and e~^ ~ 7? jCE^ , a little trial and error 
gives the following fitting function for the radiation temperature 
in terms of E and 77, 


kT, = 


Efn 

[3 - 2.449724(H'VCE3)] ’ 


(16) 


The constants in the denominator have been chosen to give the cor¬ 
rect result in the Planck and Wien limits. There is a small devia¬ 
tion at intermediate values of the chemical potential, but the frac¬ 
tional error in the temperature is at worst no more than 0.04%. In 
the present photon-conserving Comptonization approach, instead 
of using equation (|^ which corresponds to blackbody Comptoniza¬ 
tion, we use equation t o calculate T^. 

Note that the ratio E/77 is equal to 2.1012 kT^ in the black¬ 
body limit and equal to 3 kT^ in the extreme Wien limit. This is 
a small enough difference that, in light of our other approxima¬ 
tions, we could simply use the blackbody value of the ratio for all 
regimes. The key point is that the radiation temperature kT^ should 
be estimated via the ratio Efn (with some reasonable coefficient) 
if we wish to satisfy photon conservation, whereas any approach 
based purely on equation (|^ will violate photon conservation. In 
this work, we use the fitting function (T§ when we model photon- 
conserving ComptonizaitonQ 

Since we have introduced a new variable 77, we need an evolu¬ 
tion equation for this quantity. Here we use the fact that absorption 
removes photons, emission adds photons, and scattering leaves the 
number of photons unchanged. Thus, in the fluid frame, we write 
the rate of change of photon number density as 


Y k^epe 1 

'4nK^BIkTg\ 

.l3-2.449724(n‘*/C£3)j 

^ 2.701178 j 


The evolution equation for the photon number density n is then 
written as 

(nu^^).^=n, (18) 

where is the four-velocity of the radiation frame, and we have 
used the fact that 77 is a frame-invariant quantity (because h = 
dNIdt dxidx 2 dx 2 ,, and both the numerator and the denominator are 
relativistic invariants). 

In summary, under the photon-conserving Comptonization ap- 

f )ximation, we solve equations and ([TsTl using equations 

, ll^ and ([it} , where B is given by equaticmQ and is given 
by equation ( |16| ). 


3 NUMERICAL METHODS 

To solve equations |2|5| we use the general relativistic radiation 
magnetohydrodynamical (GRRMHD) code KORAL ( [S^dowski et| 


^ Because of the negative sign in the denominator of equation jl6|, the 

function can sometimes diverge and behave unphysically. This potential 
pathology (which h as n ot been an issue in any of the simulations we have 
run so far, including those presented here) can be fixed by resetting the de¬ 
nominator to its minimum value of 2.7012 whenever it tries to go lower, 
or by choosing a different form of the fitting function (Jonathan McKinney, 
private communication). 
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|al.|2013||2014| ), which employs a Godunov scheme to evolve the 
conservation equations in a fixed, arbitrary spacetime using finite- 
difference methods. The magnetic field is evolved via the induction 
equation, with the divergence-free criterion being enforced using 
a flux-constrained ( |T6th|200^ scheme as described in |Gammie et| 
[^ < |2003| ). The radiation field is evolved through its energy density 
and flux, and the radiation stress-energy tensor is closed by means 
of the Ml closure scheme ( Levermore|1984| |. For details of the ba¬ 
sic method see |S^dowski et al.| ( |2014| ). 

For the purpose of this work we adopt a simple prescription 
for the absorption opacity Ka ( |Kato et al.|2008| ), 

= 6.4 X cm^g"^ (19) 


Note that the absorption opacity should in principle be a function 
of both gas and radiation temperatures. We ignore this complication 
here. Similarly, in the spirit of keeping things simple, we ignore the 
fact that the rate of change of photon number by absorption and 
emission is generally more complicated than the prescription given 
in equation |T7] ). These issues will be dealt with more accurately in 
future studies. For the scattering opacity we use 

/^es = 0.34 cm^g"^ (20) 


Thermal Comptonization within the “blackbody” approxima¬ 
tion (Section[2^ has already been included in KORAL and has been 
described in a previous paper ( [S^dowski et al.|2015a| l. Implement¬ 
ing that scheme does not require any modification to the basic algo¬ 
rithm in KORAL; we merely need to modify the time component of 
the fluid frame four-force, which describes the rate of exchange of 
energy between gas and radiation, according to equation ( p^ . The 
corresponding radiation four-force density then takes the form 
( [S^dowski et al.|20l4) , 

G^comp. = ^Co.pc«". (21) 

where. 


'^Compt 




4k(T, - T,) 


1 -h 3.683 


nigC^ 

rriaC^ 


+ 4 


kT^ 


1 -h 


kT^ 


( 22 ) 


The photon-conserving Comptonization scheme introduced in 
this work (Section [T3] ), on the other hand, requires a modified algo¬ 
rithm. One more conserved quantity - the photon number density 
n - has to be evolved, for which we use equation ( |18| ). The source 
term on the right hand side of ( p^ describes the change of photon 
number due to emission and absorption by the gas, and is propor¬ 
tional to the absorption opacity k^. As in the case of other radiative 
source terms (equations. |4|5| ), this term may become stiff in an op¬ 
tically thick medium. Therefore, it is necessary to apply it in an 
implicit way. For this purpose we modify the original semi-implicit 
operator ( [S^dowski et al.|20T^ , which applies the standard radia¬ 
tive source terms, to include the evolution of the photon number 
density. The set of equations solved in the implicit step is now. 


(pwOo+i) - (p«')(o 

= 0, 

(23) 

'T't 'T't 

^ v,(/+l) ^ v,(i) 

= AtGy,/+l, 

(24) 

pt pt 

^V,(/+l) ^V,(/) 

= -AtGv,/+i, 

(25) 

(««;)(,+1) - 


(26) 


where subscripts (/) and (i + 1) indicate values at the beginning 
and end of the current time step, respectively. Due to the symmetry 
of the problem, it is enough to iterate just one set of primitives. 


Table 1. Model parameters 


AIbh 

lOMo 


0.0 

Pmax 

4.3 X 10“^ g/cm^ 

/^max 

10.0 

NrxNq 

252 X 234 

Rmin / ^max / Rq / Hq 

1.85/1000/0/0.6 

tmax 

48,000 

<M>/MEdd 

10.2 (blackbody Compt.) 

11.2 (N-Compt.) 

9.2 (no Compt.) 


Pmax - maximal density of the initial torus 

Pmax - maximal value of initial total to magnetic pressure ratio 

tmax - duration of simulation 


either radiative or hydrodynamical, extended by the photon number 
density, and to use the conservation of energy and momentum to 
recover the remaining quantities. 

The new system of equations is solved for values at the new 
moment of time (denoted by (/ +1)) using the Newton-Raphson 
method in five dimensions. The computational cost increases be¬ 
cause of the extra dimension by a factor of < 2 in most cases. 


4 SIMULATIONS 


4.1 Initial setup 


To assess how important the proper treatment of Comptonization 
is for accretion disks we carried out a set of simulations of super¬ 
critical accretion flows on a lOM© BH. In order to obtain a rel¬ 
atively large radial range of infiow/outfiow equilibrium, we de¬ 
cided to perfom the simulations in 2D assuming axisymmetry. The 
poloidal magnetic field was maintained by means of the mean-field 
magnetic dynamo model of |S^dowski et al.| ( [2015a| ). 

We performed three simulations, each initiated in exactly the 
same way as described in |S^dowski et al.| ( |2014| ), with parameters 
given in Table[^ The only difference between the runs was the way 
Comptonization was treated: (i) basic “blackbody Comptoniza¬ 
tion” prescription (Sectio n |2.2| ), (ii) improved “photon-conserving 
Comptonization” (Section (iii) no Comptonization at all. For 
the run with photon-conserving Comptonization we had to spec¬ 
ify the initial photon density. Because the radiation and gas in the 
initial torus satisfy local thermal equilibrium, we set the radiation 
and gas temperatures to be equal and set the photon number density 
according to the blackbody formula. 


2.7012ytB ■ 


(27) 


Hereafter, we use the gravitational radius Rg = GM/c^ as the unit 
of length, and as the unit of time. 


4.2 Accretion flow properties 

All three simulations were run up to a final time of t^ax = 48,000. 
The region of infiow/outfiow equilibrium extended up to a radius 
R ^ 30 in the equatorial plane, and to a much larger distance in 
the funnel region near the axis (where the radial velocity of the out¬ 
flowing gas is much higher, > 0.1c). The mean accretion rate for all 
the models was ~ lOMndd (Table [^. The time history of the mass 
accretion rate through the BH horizon is shown in Fig.[^ Although 
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Figure 1. Mass accretion rate at the BH as a function of time for the three Figure 3. Radial profiles of density weighted gas (solid) and radiation 

models considered. (dashed lines) temperatures in the three models. The shaded regions reflect 

the difference between the two temperatures. 


the profiles are highly variable, as expected for a turbulent accre¬ 
tion flow, there are no significant differences in the three profiles. 
This suggests that the precise treatment of Comptonization has lit¬ 
tle impact on the the gross dynamical properties of the gas. 

A similar conclusion is reached when comparing the time- and 
0-averaged distributions of density in the poloidal plane, shown 
in the top row of panels in Fig. From the left, the three pan¬ 
els correspond to blackbody Comptonization, photon-conserving 
Comptonization, and no Comptonization, respectively. In all three 
simulations, the accretion flow is geometrically thick, with an aver¬ 
age density scale height H/R ~ 0.3. The densities at the equatorial 
plane are similar in all three runs, with p ~ 10“^ g The dashed 
lines in the panels show the location of the scattering photospheres 
as measured from infinity along fixed polar angle 6. In all mod¬ 
els the polar region is optically thin (i.e., the observer looking into 
the funnel can directly see radiation coming from the vicinity of 
the BH horizon). The size of this region is largest for the photon- 
conserving model, but the difference is not large. 

Although the densities in the three runs are similar, the gas and 
radiation temperatures, shown respectively in the second and third 
rows of panels in Fig. are significantly different. In scattering 
dominated media, and in stellar mass BH accretion disks in partic¬ 
ular, Comptonization often dominates the energy transfer rate be¬ 
tween the gas and radiation; without Comptonization, the coupling 
is very weak. 

Turbulent dissipation in the disk heats the gas, the rate of heat¬ 
ing depending mostly on the properties and evolution of the MHD 
turbulence, which is insensitive to the way we treat Comptoniza¬ 
tion. When the gas subsequently cools, the dissipated heat energy 
is transfered to radiation. To sustain a given rate of energy transfer 
from gas to radiation, one can either have efficient coupling, with 
gas and radiation temperatures close to each other, or weak cou¬ 
pling, with a large difference in the two temperatures. In the case 
of the simulation with no Comptonization (right panels), the cou¬ 
pling is weak, therefore the gas temperature becomes very large. 
This effect is significant in the inner regions of the disk and is par¬ 
ticularly evident in the coronal and polar regions, where the gas 
temperature exceeds 10^^ K, while the radiation temperature stays 
around 10^ K. In contrast, for the two models that allow Compton 
coupling, viz., blackbody Comptonization and photon-conserving 
Comptonization, the gas temperature is much lower, and is also 


closer to the radiation temperature, reflecting the fact that the cou¬ 
pling between the gas and the radiation field is strong. As in the 
case of the no Comptonization model, the gas temperature tends 
to be higher than the radiation temperature here as well, but the 
difference, even in the coronal region, is no longer very extreme. 

Radial profiles of the gas and radiation temperatures (these 
are density-weighted, and therefore are dominated by the gas in the 
disk interior rather than the corona) are shown in Fig.[^ The aver¬ 
age gas temperature at radius R - 15 is approximately 5 x 10^, 
6 X 10^, and 3 x 10^ K for models with no, photon-conserving, 
and blackbody Comptonization, respectively. The overheating of 
the gas in the no Comptonization model is dramatically evident. 
The other two models (which include Comptonization) have lower 
gas and radiation temperatures and the two temperatures are nearly 
equal in the bulk of the disk {R> 10). The actual temperatures are, 
however, not the same for the two models. The photon-conserving 
model is hotter by a factor of 2 (measured by either gas or radiation 
temperature). This is a significant difference and will have a notice¬ 
able effect on the spectral properties of the escaping radiation. 

The blackbody Comptonization and no Comptonization mod¬ 
els assume that the radiation spectrum is blackbody. Therefore, in 
these models the radiation temperature is a function only of the 
radiative energy density in the fluid frame. 



For the photon-conserving Comptonization model, however, we 
track the number density of photons and use equation |T6t to com¬ 
pute T^. This temperature is equal to Tj -bb only in those regions of 
the flow where the radiation happens to be black body, e.g., deep 
inside the disk. In other regions, the local number density of pho¬ 
tons deviates from the blackbody value, and so do the temperature 
and spectrum. We define the color correction factor as the ratio of 
the real radiation temperature to the blackbody temperature. 


Values of fco\ larger than unity correspond to a radiation spectrum 
that is harder than blackbody. Under the gray approximation we 
have adopted in this work, we do not have detailed information on 
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Figure 2. Distributions of density (top row), gas temperature (second row) radiation temperature (third row), and cooling rate (bottommost row) for models 
with blackbody Comptonization (left), photon-conserving Comptonization (middle), and no Comptonization (right panels). Plots correspond to data averaged 
over t = 24000 ^ 48000. The dashed lines in the density plots reflect the location of the scattering photosphere as measured from inflnity along flxed polar 
angle. The contours on the radiation temperature plots (third row) reflect the radiation to gas temperature ratio. 


(c) 0000 RAS, MNRAS 000, 000-000 













































Comptonization in accretion disks 1 



0 10 20 30 40 50 


rjU 

Figure 4. Color correction factor /coi (equation |29) corresponding to the 
photon-conserving Comptonization model. The colors reflect the distribu¬ 
tion of gas density and the black dashed line shows the location of the ef¬ 
fective photosphere (equation!^. 


the spectral shape of the radiation field. Therefore, spectral distor¬ 
tion from blackbody is described by a single number /coi. More 
detailed, frequency-dependent radiative transfer is required if we 
wish to obtain additional information. 

By definition, the color correction factor (equation]^ is equal 
to unity for the blackbody and no-Comptonization models, since 
the radiation spectrum is assumed to be blackbody. In the photon- 
conserving Comptonization model, however, it can (and does) de¬ 
viate from unity. The contours in Fig. show the color correction 
factor for this model plotted on top of the gas density distribution. 
The black dashed line shows the location of the effective photo¬ 
sphere, estimated by integrating the effective optical depth from 
the axis along Q at fixed radius and identifying the angle at which 
it equals 2/3: 

P + 'fes) dff = 113. (30) 

Deep in the disk and at large radii {R > 30), where the disk is 
opticaly thick to both scattering and absorption, f^oi is close to unity 
and the radiation spectrum is close to blackbody. The color correc¬ 
tion increases (the spectrum becomes harder) as the radiation gets 
closer to the disk surface — the photons are up-scattered by hot 
gas and gain energy. Near the edge of the funnel, the color correc¬ 
tion factor reaches f^oi ~ 5.0, which reflects how strongly scatter¬ 
ing modifies the radiation fieldj^This large value suggests that the 
spectrum would be strongly dominated by a Compton hump, and 
the original black body component would be insignificant (com¬ 
pare |Kawashima et al.|2012| l . 

The rate of energy transfer from gas to radiation is the cooling 
rate of the gas and describes how effectively the heat energy re¬ 
leased by turbulent dissipation goes into the radiation field. The 
bottommost panels of Fig. show the distribution of the time- 
averaged cooling rate for the three models. For all three simula¬ 




r/M 

Figure 5. Radiative luminosity Fj-ad (equation!^ as a function of radius 
for the three models. 


tions, the highest rates of energy transfer are at the equatorial plane 
and in the innermost region (R < 15). However, the vertical extent 
of this region as well as the magnitude of the cooling differ. The no 
Comptonization model shows the least amount of cooling both in 
the funnel region near the axis and at the equatorial plane. For the 
thermal and photon-conserving Comptonization models, the cool¬ 
ing region extends further towards the axis and the magnitude of 
the cooling rate is also larger. 

Because accretion flows are quasi-stationary, the energy trans- 
fered from gas to radiation via cooling must be removed in some 
fashion. In optically thick disks, this can happen in two ways — 
either by radiation flowing out of the disk or by radiation being 
trapped in the gas and carried into the BH with the flow (advected). 
Fig.j^is a plot of the radiative luminosity as a function of radius for 
the three models, calculated as follows, 

£rad = J R^,sFgdA, (31) 

where is the metric determinant and dA = 2nr^d6 is a dif¬ 
ferential surface element at radius r. Negative values of this quan¬ 
tity correspond to regions where radiation advection into the BH 
dominates. The radius where £rad = 0 reflects the location of the 
effective trapping radius of the accretion flow. 

All three runs show significant advection of photons. Only 
outside radius R ^ 25 does the net flux of radiation point outward. 
To estimate the total amount of energy put into the radiation field 
we take the luminosity in inflowing photons at the BH and add it 
to the luminosity in outflowing photons at radius = 45. The to¬ 
tal energy release in radiation is roughly 7%, 6%, and 4%Mc^ for 
the blackbody, photon-conserving, and no Comptonization models, 
respectivel}|^ These numbers are in agreement with the discussion 
of the cooling rate in the previous paragraphs, viz., not including 
Comptonization underestimates the cooling rate. 

By taking the radial derivative of the radiation luminosity 
we obtain the local photon generation rate or the gas cooling rate. 
This quantity is shown in Fig.|^ The emission profiles more or less 
overlap in the outer (R > 15) region but are significantly different 


^ It is interesting to note that similar values of the color correction factor 
were obtained by |Kawaguchi| {200^ who estimated Comptonized spectra 
of slim accretion disks. 


^ We stress here that these numbers do not correspond to the standard def¬ 
inition of radiative efficiency, which includes only radiation escaping to a 
distant observer. 
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Figure 6. Profiles of radiative emission per unit radius for the three models. 

inside that radius. The maximum amount of energy is put into radi¬ 
ation in the blackbody Comptonization model, and the least when 
Comptonization is not taken into account. In the latter model, there 
is hardly any emission inside the marginally stable orbit. 


5 DISCUSSION 

In this work we have introduced an improved method for treating 
Comptonization in computational fluid dynamics. The suggested 
approach conserves the number of photons and allows for a non- 
blackbody spectrum of radiation. It is fully covariant and may be 
applied to simulations performed in general relativity. However, the 
method is still very approximate, since it assumes a specific spectral 
shape for the radiation, viz., a Bose-Einstein distribution function 
(eq. [ID- Also, the proposed algorithm is valid only for gray ra¬ 
diative transfer. Nevertheless, we believe this “photon-conserving 
Comptonization” method is a significant improvement over the tra¬ 
ditional “blackbody Comptonization method”. 

To assess the effect of various treatments of Comptonization 
in modeling accretion disks, we performed a set of three simula¬ 
tions corresponding to a super-critical disk accreting at lOMEdd- 
We implemented the new photon-conserving approach (described 
in detail in Section [0] ) as well as blackbody Comptonization (Sec¬ 
tion [2^ and considered also a model with no Comptonization. We 
And that the density distribution and the accretion rate are insen¬ 
sitive to the way we model Comptonization, whereas the gas tem¬ 
perature and the properties of the radiation held depend strongly 
on which version of Comptonization we use. In particular, simula¬ 
tions with no Comptonization produce unphysically high gas tem¬ 
peratures. This is because of inefficient coupling between gas and 
radiation. The blackbody approach, on the other hand, noticeably 
underestimates the gas and radiation temperatures. In making these 
statements, we are working on the assumption that the photon- 
conserving model, being the most sophisticated of the three, is clos¬ 
est to the true solution. 

By comparing the three simulations we find that the impact of 
Comptonization is not limited to the coronal regions (as claimed 
by some authors), but that it actually affects the gas properties 
throughout the disk, even deep inside the effective photosphere. 
With no Comptonization at all, the coupling between gas and ra¬ 
diation comes only from the absorption opacity, which requires 


significantly different gas and radiation temperatures. With Comp¬ 
tonization, the optically thick gas in the bulk of the disk is at lower 
temperatures, closer to the radiation temperature. 

The photon-conserving model provides one additional piece 
of information compared to the other two models, viz., the color 
correction factor f^oi as a function of position. This factor, which 
measures the ratio of the radiation temperature to an effective 
blackbody temperature (eq. [^, describes the amount of spectral 
hardening caused by scattering and Comptonization. We find that 
/col can deviate substantially from unity, and can be as large as 5 
or 6 at the photosphere inside the funnel. Such high color correc¬ 
tion factors will characterize radiation escaping through the polar 
funnel region, and will be visible to observers viewing at small po¬ 
lar angles. For larger inclinations, most of the observed radiation 
would come from the disk photosphere located far from the BH 
(this region is not resolved in our simulations) and would likely not 
show such strong Comptonization. 

Another quantity that is significantly affected by the version 
of Comptonization that one uses is the amount of cooling, which 
determines the ultimate luminosity (both in advected and escaping 
photons). The blackbody and no-Comptonization models over- and 
under-estimate the cooling rate, respectively. 

Taking into account the modest computational cost of the 
photon-conserving Comptonization approach, we recommend us¬ 
ing it in simulations of accretion disks. It is true that blackbody 
Comptonization reproduces reasonably well the dynamical proper¬ 
ties of the accretion flow and the cooling rate. However, the gas 
and radiation temperatures obtained via this approach are signifi¬ 
cantly under-estimated. Therefore, this model cannot be trusted for 
applications in which the temperature is of interest, e.g., spectral 
modeling. 

To complete this study we need to assess how the treatment 
of Comptonization affects the spectrum as seen by a distant ob¬ 
server. This requires post-processing with a multi-dimensional ra¬ 
diative transfer solver that works in general relativity and allows 
for frequency-dependent opacity (e.g., |Zhu et al.|20I5] l and Comp¬ 
tonization. 
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